{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../common/rfsoc_book_banner.jpg\" alt=\"University of Strathclyde\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook Set D\n",
    "\n",
    "---\n",
    "\n",
    "## 03 - Pulse Shaping\n",
    "In this notebook, we will cover some basic pulse shaping techniques such as square pulse shaping, sinc and raised cosine pulse shapes, matched filtering and the minimization of Inter Symbol Interference (ISI).\n",
    "\n",
    "## Table of Contents\n",
    "* [1. Introduction](#introduction)\n",
    "* [2. Transmitting Pulses](#pulses)\n",
    "* [3. Square Pulse Shaping](#squares)\n",
    "    * [3.1. Averaging at the Receiver](#averaging-at-the-receiver)\n",
    "* [4. Sinc Function](#sinc)\n",
    "* [5. Raised Cosine](#rcos)\n",
    "    * [5.1. Raise Cosine Pulse Shaping](#raised-cosine-pulse-shaping)\n",
    "* [6. Matched Filtering with Root Raised Cosine](#rrc)\n",
    "* [7. Conclusion](#conclusion)\n",
    "\n",
    "## Revision\n",
    "* **v1.0** | 05/12/22 | *First Revision*\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction <a class=\"anchor\" id=\"introduction\"></a>\n",
    "\n",
    "Pulse shaping is an integral part of many communication systems. It serves as a means to band-limit the baseband transmit signal and to improve robustness to noise. Usually, pulse shaping is accomplished by passing our symbols (pulses) through a pulse shaping filter.\n",
    "\n",
    "You will only need the Numpy and Matplotlib libraries to run this notebook. Begin by importing the necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Transmitting Pulses <a class=\"anchor\" id=\"pulses\"></a>\n",
    "Previously in the baseband modulation notebook, we simulated BPSK (Binary Phase Shift Keying) transmissions, where each symbol was a pulse and represented a bit. However, if we try to transmit a pulse like that over the air we would cause a lot of problems to surrounding radio devices through interference - a pulse (or in this case an impulse) in the time domain actually contains all frequencies in the frequency domain.\n",
    "\n",
    "We will use the functions defined in previous notebooks to generate some BPSK symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to add white gaussian noise to a signal, specify SNR in dB\n",
    "def awgn(signal, snr):\n",
    "    sig_power = np.mean(np.abs(signal)**2) # calculate signal power\n",
    "    sig_power_db = 10* np.log10(sig_power) # convert to dB\n",
    "\n",
    "    noise_power_db = sig_power_db - snr\n",
    "    noise_power = 10**(noise_power_db / 10)\n",
    "\n",
    "    complex_noise = np.sqrt(noise_power/2)*(np.random.randn(len(signal)) + \\\n",
    "                                            np.random.randn(len(signal))*1j)\n",
    "\n",
    "    return signal + complex_noise\n",
    "\n",
    "# Function to generate BPSK\n",
    "def generate_bpsk(num_symbols, noise=50):\n",
    "    bits = np.random.randint(0,2,num_symbols)\n",
    "    bpsk_scheme = [1+0j, -1+0j]\n",
    "    bpsk_symbols = np.array([bpsk_scheme[i] for i in bits])\n",
    "    \n",
    "    bpsk_symbols = awgn(bpsk_symbols, noise)\n",
    "    \n",
    "    return bpsk_symbols"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now generate BPSK symbols and plot them in the time domain. BPSK only uses the real channel for binary signalling, so we disregard the imaginary part."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "symbols = generate_bpsk(2048)\n",
    "plt.stem(symbols.real[:8], use_line_collection=True)\n",
    "plt.xlabel('Samples')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.title('BPSK Pulses')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While this may seem a straightforward way to transfer bits around, looking at the frequency domain of these pulses we can see that they contain all frequencies. This would be undesirable, because it would cause interference to other channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 128\n",
    "fft_size = 2048\n",
    "\n",
    "y_fft = np.fft.fftshift(np.fft.fft(generate_bpsk(fft_size),fft_size))[int(fft_size/2):]\n",
    "freqs = np.fft.fftshift(np.fft.fftfreq(fft_size,1/fs))[int(fft_size/2):]\n",
    "\n",
    "plt.plot(freqs, np.abs(y_fft))\n",
    "plt.xlabel('Frequency, Hz')\n",
    "plt.ylabel('Magnitude')\n",
    "plt.title('BPSK Pulses in Frequency Domain')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Square Pulse Shaping <a class=\"anchor\" id=\"squares\"></a>\n",
    "\n",
    "One of the simplest forms of pulse shaping is the square pulse shape, where we repeat the same pulse $N$ times per symbol - usually we refer to $N$ as samples per symbol, or sps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of symbol repeats\n",
    "N = 8\n",
    "\n",
    "# Apply square shape to the BPSK symbols\n",
    "symbols_square_shaped = np.repeat(symbols,N)\n",
    "\n",
    "# Plot some of the pulse shaped symbols\n",
    "plt.stem(symbols_square_shaped.real[:64], use_line_collection=True)\n",
    "plt.xlabel('Samples')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While this is better than transmitting pulses, square shapes still produce a lot of spectral emissions due to the sharp transitions in the time domain. If we take a clean square wave as an example, we can see that it is made of a combination of many harmonics in the frequency domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create clean square wave\n",
    "square_wave = np.repeat(np.array([0, 1, 0, 1, 0, 1, 0, 1, 0, 1]),8)\n",
    "\n",
    "# Take the fft of the square wave and calculate corresponding x-axis frequencies\n",
    "y_fft = np.fft.fftshift(np.fft.fft(square_wave,fft_size))[int(fft_size/2):]\n",
    "freqs = np.fft.fftshift(np.fft.fftfreq(fft_size,1/fs))[int(fft_size/2):]\n",
    "\n",
    "fig, ax = plt.subplots(1,2, figsize=(10,4))\n",
    "ax[0].plot(square_wave)\n",
    "ax[0].set_xlabel('Samples')\n",
    "ax[0].set_ylabel('Amplitude')\n",
    "ax[0].set_title('Time Domain')\n",
    "ax[1].plot(freqs, np.abs(y_fft))\n",
    "ax[1].set_xlabel('Frequency, Hz')\n",
    "ax[1].set_title('Frequency Domain')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a similar way, the FFT of a square pulse shaped waveform will create many side lobes. These have the potential to interfere with users on adjacent frequencies, once the signal is modulated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_fft = np.fft.fftshift(np.fft.fft(symbols_square_shaped,fft_size))[int(fft_size/2):]\n",
    "freqs = np.fft.fftshift(np.fft.fftfreq(fft_size,1/fs))[int(fft_size/2):]\n",
    "\n",
    "plt.plot(freqs, np.abs(y_fft))\n",
    "plt.xlabel('Frequency, Hz')\n",
    "plt.ylabel('Magnitude')\n",
    "plt.title('Square Shaped Pulses in Frequency Domain')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Averaging at the receiver <a class=\"anchor\" id=\"averaging-at-the-receiver\"></a>\n",
    "A square pulse shape receiver can be implemented using a moving average filter. We use 4 weights and set them to [0.25, 0.25, 0.25, 0.25], as we are oversampling with $N=4$ samples per symbol. We anticipate the filter to output a peak response of 1 once the filter fully overlaps with an incoming square pulse shape of 4 samples. Note that using a square pulse shape filter at the transmitter and a moving average filter at the receiver, we do achieve a form of matched filtering, because both filters share the same square shaped response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = [0.25, 0.25, 0.25, 0.25]\n",
    "filtered_symbols = np.convolve(np.repeat(symbols,4).real, weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convolving this filter across our pulse shaped signal will create triangular peaks or 'maximum effect points' that will minimise the chance of incorrect symbol decisions (under noisy conditions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.stem(filtered_symbols.real[:64], use_line_collection=True)\n",
    "plt.plot(filtered_symbols.real[:64])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting our filtered signal we can see that the side lobes have been suppressed, and as a result we are experience less interference at the receiver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the one-sided FFT of the filtered signal\n",
    "y_fft = np.fft.fftshift(np.fft.fft(filtered_symbols,fft_size))[int(fft_size/2):]\n",
    "freqs = np.fft.fftshift(np.fft.fftfreq(fft_size,1/fs))[int(fft_size/2):]\n",
    "\n",
    "# Plot the result\n",
    "plt.plot(freqs, np.abs(y_fft))\n",
    "plt.xlabel('Frequency, Hz')\n",
    "plt.ylabel('Magnitude')\n",
    "plt.title('Square Shaped Pulses in Frequency Domain')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spectral properties of this receiver are better than transmitting impulses and consuming unnecessary bandwidth. Square pulse shapes still introduce a lot of interference to nearby frequencies. Surely we can do better!\n",
    "\n",
    "## 4. Sinc Function <a class=\"anchor\" id=\"sinc\"></a>\n",
    "The ideal bandlimiting response for transmitting our signal would be a “brick wall” — i.e. perfectly passing frequencies up to the actual signal bandwidth, with zero emissions on adjacent frequencies. This \"brick wall\" response can be achieved by using a sinc shaped filter in the time domain.\n",
    "\n",
    "The sinc function is given by\n",
    "\n",
    "$$\n",
    "n(t) = \\frac{sin(\\pi t / T)}{\\pi t / T}\n",
    "$$\n",
    "\n",
    "A true sinc function is infinite. However, this is not realisable, as we would need an infinitely expensive filter with an infinite number of weights. In practice, a truncated sinc shape filter is used instead. Much like windowing functions, it has outputs that taper off at the ends and smooth out the response in order to reduce the sharp transitions that cause unwanted frequency emissions.\n",
    "\n",
    "In Numpy we can simply call `np.sinc`.\n",
    "\n",
    "Lets visualize a Sinc shape filter that is 40 symbol periods long and uses 5 samples per symbol. To realize this filter, we will need 251 weights (40 symbols $\\times$ 5 samples per symbol + 1 center weight = 251). We assume a sample period of $T$=1s (sampling rate of 1Hz) - in the next code cell we will omit $T$ because it is just 1. Just keep in mind that it is part of the formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 1 # sample period\n",
    "sps = 5 # symbol period in samples\n",
    "num_weights = 251\n",
    "x = np.arange(-int(num_weights/2),int(num_weights/2)+1,1)/sps/T\n",
    "sinc_weights = np.sinc(x)\n",
    "\n",
    "plt.plot(x,sinc_weights)\n",
    "plt.title('Sinc function')\n",
    "plt.xlabel('Time samples, t')\n",
    "plt.ylabel('n(t)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before inspecting the sinc response in the frequency domain, let's generate another filter that runs for only 4 symbol periods to compare the frequency responses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create another, shorter sinc filter\n",
    "# num_weights = 4 symbol periods * 5 samples per symbol + 1\n",
    "num_weights = 21\n",
    "x = np.arange(-int(num_weights/2),int(num_weights/2)+1,1)/sps\n",
    "sinc_weights2 = np.sinc(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing the filter responses in the frequency domain, we can see that the more symbol samples the filter is allowed, the closer we get to the ideal cut-off. Desired frequencies of the main lobe are passed, while any unnecessary emissions are suppressed. This relationship basically indicates that more filter weights in time domain means better spectral properties in frequency domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate frequency responses of the 2 filters\n",
    "fft_size = 1024\n",
    "sinc1_response = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft(sinc_weights, fft_size)))/fs)\n",
    "sinc2_response = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft(sinc_weights2, fft_size)))/fs)\n",
    "\n",
    "# Plot the frequency responses\n",
    "plt.plot(sinc1_response[int(fft_size/2):-250]) # truncated for better viewing experience\n",
    "plt.plot(sinc2_response[int(fft_size/2):-250])\n",
    "\n",
    "plt.legend(('50 symbol periods', '4 symbol periods'))\n",
    "plt.title('Sinc filter magnitude responses')\n",
    "plt.ylabel('Magnitude, Hz')\n",
    "plt.xlabel('Frequency, Hz')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we are limited in the number of filter weights we want to use (which we often are), but still want to minimze the outer frequency emissions? Well we can improve our sinc filter.\n",
    "\n",
    "## 5. Raised Cosine <a class=\"anchor\" id=\"rcos\"></a>\n",
    "The Raised Cosine (RC) pulse shape has a preferable impulse response that tails off more quickly (as governed by the design parameter α). It is defined as:\n",
    "\n",
    "$$\n",
    "p(t) = (\\frac{sin(\\pi t / T)}{\\pi t / T}) (\\frac{cos ( \\alpha \\pi t / T)}{1 - (2\\alpha t / T)^{2}})\n",
    "$$\n",
    "\n",
    "Note that if $\\alpha = 0$, the second part of the equation will be equal to 1, and we end up with a response equivalent to the sinc function.\n",
    "\n",
    "Here, we can visualize what effect the parameter $\\alpha$ has on our new Raised Cosine filter response. We can see that the cosine window term suppresses the tails of the Sinc shape as the $\\alpha$ term increases - this allows us to design a cheaper pulse shaping filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sps = 5\n",
    "num_weights = 101\n",
    "x = 0.9999*np.arange(-int(num_weights/2),int(num_weights/2)+1,1)/sps\n",
    "\n",
    "alphas = [0, 0.25, 0.5, 0.75]\n",
    "for alpha in alphas:\n",
    "    raised_cos_weights = np.sinc(x)*(np.cos(alpha*np.pi*x)/(1-((2*alpha*x)**2)))\n",
    "    plt.plot(x,raised_cos_weights)\n",
    "\n",
    "plt.legend((r'$\\alpha = 0$', r'$\\alpha = 0.25$', r'$\\alpha = 0.5$', r'$\\alpha = 0.75$'))\n",
    "plt.ylabel('Filter weight magnitude')\n",
    "plt.xlabel('Filter weight index')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we look at the filter responses at different $\\alpha$ values, we can see that the required bandwidth to capture our signal increases with $\\alpha$. Like most things, the $\\alpha$ parameter is a tradeoff between how much spectral leakage you are allowing to nearby frequencies and how much bandwidth will be required to recover the transmitted signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sps = 5\n",
    "num_weights = 121\n",
    "x = 0.9999*np.arange(-int(num_weights/2),int(num_weights/2)+1,1)/sps\n",
    "\n",
    "for alpha in alphas:\n",
    "    raised_cos_weights = np.sinc(x)*(np.cos(alpha*np.pi*x)/(1-((2*alpha*x)**2)))\n",
    "\n",
    "    rcos_response = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft(raised_cos_weights, fft_size)))/sps)\n",
    "    plt.plot(rcos_response[int(fft_size/2):-100])\n",
    "    \n",
    "plt.legend((r'$\\alpha = 0$', r'$\\alpha = 0.25$', r'$\\alpha = 0.5$', r'$\\alpha = 0.75$'))\n",
    "plt.title('Magnitude Response')\n",
    "plt.xlabel('Frequency, Hz')\n",
    "plt.ylabel('Magnitude, dB')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1. Raised Cosine Pulse Shaping <a class=\"anchor\" id=\"raised-cosine-pulse-shaping\"></a>\n",
    "\n",
    "We will begin by using a simple 8 symbol period, or 41 weight RC filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 5\n",
    "num_weights = 41\n",
    "alpha = 0.5\n",
    "x = 0.9999*np.arange(-int(num_weights/2),int(num_weights/2)+1,1)/fs\n",
    "raised_cos_weights = np.sinc(x)*(np.cos(alpha*np.pi*x)/(1-((2*alpha*x)**2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we generate 3 pulses, each 5 samples or 1 symbol period apart, so they do not interfere with each other. Notice how pulse shaping causes a delay in the signal - our first peak is no longer at index 5, but 25! This is because pulse shaping introduces latency, and for our RC filter this is half the total pulse duration, which is 4 symbol periods (or 20 samples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate pulse vectors\n",
    "pulse1 = np.zeros(21,)\n",
    "pulse2 = np.zeros(21,)\n",
    "pulse3 = np.zeros(21,)\n",
    "pulse1[5] = 1\n",
    "pulse2[10] = 1\n",
    "pulse3[15] = 1\n",
    "\n",
    "# Pulse shape each symbol separately \n",
    "symbol1_rc = np.convolve(pulse1,raised_cos_weights).real\n",
    "symbol2_rc = np.convolve(pulse2,raised_cos_weights).real\n",
    "symbol3_rc = np.convolve(pulse3,raised_cos_weights).real\n",
    "\n",
    "# Pulse shape the combination of symbols\n",
    "symbols_rc = np.convolve(pulse1+pulse2+pulse3,raised_cos_weights).real\n",
    "\n",
    "# Plot them all\n",
    "fig, ax = plt.subplots(1,3, figsize=(16,4))\n",
    "\n",
    "ax[0].plot(pulse1, '-o')\n",
    "ax[0].plot(pulse2, '-o')\n",
    "ax[0].plot(pulse3, '-o')\n",
    "ax[0].set_title('Individual Pulses')\n",
    "ax[0].set_xlabel('Samples')\n",
    "ax[0].set_ylabel('Amplitude')\n",
    "ax[0].legend(('Pulse 1', 'Pulse 2', 'Pulse 3'), loc='center left')\n",
    "\n",
    "ax[1].plot(symbol1_rc)\n",
    "ax[1].plot(symbol2_rc)\n",
    "ax[1].plot(symbol3_rc)\n",
    "ax[1].plot([15,20,25,30,35,40,45], [0,0,0,0,0,0,0], 'x')\n",
    "ax[1].set_title('RC Pulse Shaped Symbols')\n",
    "ax[1].set_xlabel('Samples')\n",
    "ax[1].set_ylabel('Amplitude')\n",
    "ax[1].legend(('Symbol 1', 'Symbol 2', 'Symbol 3', 'Zero Crossings'), loc='center left')\n",
    "\n",
    "ax[2].plot(symbols_rc)\n",
    "ax[2].plot([5+20,10+20,15+20], [1,1,1],'x')\n",
    "ax[2].set_title('RC Pulse Shaped Transmission')\n",
    "ax[2].set_xlabel('Samples')\n",
    "ax[2].set_ylabel('Amplitude')\n",
    "ax[2].legend(('Transmission', 'Maximum Effect Points'), loc='center left')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The great thing about the RC filter is that it does not introduce Inter Symbol Interference (ISI) - that is the maximum effect points are unaffected by the neighbouring symbols. While they do overlap with each other, their amplitude is exactly 0 at the peaks of other symbols.\n",
    "\n",
    "## 6. Matched Filtering with Root Raised Cosine <a class=\"anchor\" id=\"rrc\"></a>\n",
    "\n",
    "In order to reap the benefits of the raised cosine with matched filtering, we split it into two filters with identical responses. The RC shape is split into two, resulting in two (identical, or “matched”) Root Raised Cosine (RRC) filters. The RRC filter is one of the most popular matched filtering options in many communications systems. The mathematics of taking a square root in the frequency domain are fairly complicated. It is not a trivial filter to implement in the time domain, as can be seen by the impulse response definition in the equation below.\n",
    "\n",
    "$$\n",
    "h(t) =\n",
    "    \\begin{cases}\n",
    "      \\frac{1}{T_s}(1+\\alpha(\\frac{4}{\\pi}-1)), & \\text{if } t=0,\\\\\n",
    "      \\frac{\\alpha}{T_{s}\\sqrt{2}}[(1+\\frac{2}{\\pi})sin(\\frac{\\pi}{4\\alpha}) + (1 - \\frac{2}{\\pi})cos(\\frac{\\pi}{4\\alpha})] & \\text{if } t=\\pm\\frac{T_s}{4\\alpha},\\\\\n",
    "      \\frac{1}{T_s}\\frac{sin[\\pi\\frac{t}{Ts}(1-\\alpha)]+4\\alpha\\frac{t}{T_s}cos[\\pi\\frac{t}{T_s}(1+\\alpha)]}{\\pi\\frac{t}{T_s}[1-(4\\alpha\\frac{t}{T_s})^2]} & \\text{otherwise.}\n",
    "    \\end{cases}\n",
    "$$\n",
    "    \n",
    "Luckily most filter design tools come equipped with RRC implementations, due to it being such a popular choice for matched filtering. If you were using the filter design tools like those found in GNU Radio or MATLAB, you need not worry about these types of equations. Here, we will implement a simplified version, by making a small filter following the 3rd part of this equation, which should be enough to demonstrate the effectiveness of this type of filtering. We will see a division error, because we are applying our formula to a value that is undefined. We compensate for this by defining our middle filter weight as shown in the 1st part of the equation. The term $\\frac{t}{T_s}$ in our code cell is defined as $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n",
    "\n",
    "# Base parameters, same as rcos in the previous example\n",
    "fs = 5\n",
    "Ts =  1/fs\n",
    "alpha = 0.5\n",
    "\n",
    "num_weights = 41\n",
    "x = np.arange(-int(num_weights/2),int(num_weights/2)+1,1)/fs\n",
    "\n",
    "# Impulse response according to definition\n",
    "h_rrc = 1/Ts*(np.sin(np.pi*x*(1-alpha)) + 4*alpha*x*np.cos(np.pi*x*(1+alpha)))/(np.pi*x*(1-(4*alpha*x)**2))\n",
    "\n",
    "# Find and replace the center weight according to the first part of the formula\n",
    "h_rrc[int(num_weights/2)] = 1/Ts*(1+alpha*(4/np.pi - 1))\n",
    "\n",
    "# Normalize the weights\n",
    "h_rrc = h_rrc/np.max(h_rrc)\n",
    "\n",
    "# Plot filter impulse response\n",
    "plt.plot(h_rrc)\n",
    "plt.title('Raised Root Cosine Impulse Response')\n",
    "plt.xlabel('Time samples, t')\n",
    "plt.ylabel('Filter response, h(t)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike the raised cosine, the RRC does introduce ISI to our transmitted pulses, which is why it is so important to pair it with the same filter on the receiver. If we pass our pulses through this filter a second time, the zero-ISI property is restored and we can see that the maximum effect points match the zero crossings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate pulse vectors\n",
    "pulse1 = np.zeros(21,)\n",
    "pulse2 = np.zeros(21,)\n",
    "pulse3 = np.zeros(21,)\n",
    "pulse1[5] = 1\n",
    "pulse2[10] = 1\n",
    "pulse3[15] = 1\n",
    "\n",
    "# Pulse shape each symbol separately \n",
    "symbol1_rrc = np.convolve(pulse1,h_rrc).real\n",
    "symbol2_rrc = np.convolve(pulse2,h_rrc).real\n",
    "symbol3_rrc = np.convolve(pulse3,h_rrc).real\n",
    "\n",
    "# Pulse shape with matched filter\n",
    "symbol1_matched = np.convolve(symbol1_rrc,h_rrc).real[20:-20]\n",
    "symbol2_matched = np.convolve(symbol2_rrc,h_rrc).real[20:-20]\n",
    "symbol3_matched = np.convolve(symbol3_rrc,h_rrc).real[20:-20]\n",
    "\n",
    "# Plot them all\n",
    "fig, ax = plt.subplots(1,2, figsize=(12,4))\n",
    "\n",
    "# Only RRC filtered pulses\n",
    "ax[0].plot(symbol1_rrc)\n",
    "ax[0].plot(symbol2_rrc)\n",
    "ax[0].plot(symbol3_rrc)\n",
    "ax[0].plot([15,20,25,30,35,40,45], [0,0,0,0,0,0,0], 'x')\n",
    "ax[0].set_title('RRC Pulse Shaped Symbols')\n",
    "ax[0].set_ylabel('Magnitude')\n",
    "ax[0].set_xlabel('Samples')\n",
    "ax[0].legend(('Symbol 1', 'Symbol 2', 'Symbol 3', 'Zero Crossings'), loc='center left')\n",
    "\n",
    "# Twice RRC filtered pulses, demonstrating matched filtering\n",
    "ax[1].plot(symbol1_matched)\n",
    "ax[1].plot(symbol2_matched)\n",
    "ax[1].plot(symbol3_matched)\n",
    "ax[1].plot([15,20,25,30,35,40,45], [0,0,0,0,0,0,0], 'x')\n",
    "ax[1].set_title('After Matched Filtering')\n",
    "ax[1].set_ylabel('Magnitude')\n",
    "ax[1].set_xlabel('Samples')\n",
    "ax[1].legend(('Symbol 1', 'Symbol 2', 'Symbol 3', 'Zero Crossings'), loc='center left')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Conclusion\n",
    "\n",
    "In this notebook we covered some pulse shaping essentials, including the sinc and raised cosine pulse shapes. We showcased the tradeoffs of filter size and desired frequency response, as well as how the $\\alpha$ parameter for the raised cosine filter can allow us to tune some spectral properties of the pulse shaping filter without changing the number of available filter weights. The next set of notebooks will discuss amplitude modulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "[⬅️ Previous Notebook](02_evm_and_ber.ipynb) || [Next Notebook 🚀](../notebook_E/01_amplitude_modulation.ipynb)\n",
    "\n",
    "Copyright © 2023 Strathclyde Academic Media\n",
    "\n",
    "---\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
